IMPROVED LINEAR RESPONSE FOR STOCHASTICALLY 

DRIVEN SYSTEMS 



RAFAIL V. ABRAMOV 

Abstract. The recently developed short-time linear response algorithm, 
which predicts the average response of a nonlinear chaotic system with 
forcing and dissipation to small external perturbation, generally yields 
high precision of the response prediction, although suffers from numer- 
ical instability for long response times due to positive Lyapunov expo- 
nents. However, in the case of stochastically driven dynamics, one typ- 
ically resorts to the classical fluctuation-dissipation formula, which has 
the drawback of explicitly requiring the probability density of the statisti- 
cal state together with its derivative for computation, which might not be 
available with sufficient precision in the case of complex dynamics (usu- 
ally a Gaussian approximation is used). Here we adapt the short-time 
linear response formula for stochastically driven dynamics, and observe 
that, for short and moderate response times before numerical instability 
develops, it is generally superior to the classical formula with Gaussian 
approximation for both the additive and multiplicative stochastic forcing. 
Additionally, a suitable blending with classical formula for longer re- 
sponse times eliminates numerical instability and provides an improved 
response prediction even for long response times. 



1. Introduction 

The fluctuation-dissipation theorem (FDT) is one of the cornerstones of 
modern statistical physics. Roughly speaking, the fluctuation-dissipation 
theorem states that for dynamical systems at statistical equilibrium the av- 
erage response to small external perturbations can be calculated through 
the knowledge of suitable correlation functions of the unperturbed dynam- 
ical system. The fluctuation-dissipation theorem has great practical use in 
a variety of settings involving statistical equilibrium of baths of identical 
gas or liquid molecules, Ornstein-Uhlenbeck Brownian motion, motion of 
electric charges, turbulence, quantum field theory, chemical physics, phys- 
ical chemistry and other areas. The general advantage provided by the 
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fluctuation-dissipation theorem is that one can successfully predict the re- 
sponse of a dynamical system at statistical equilibrium to an arbitrary small 
external perturbation without ever observing the behavior of the perturbed 
system, which offers great versatility and insight in understanding behav- 
ior of dynamical processes near equilibrium in numerous scientific applica- 
tions |fTUl[T6ll . In particular, there has been a profound interest among the 
atmospheric/ocean science community to apply the fluctuation-dissipation 
theorem to predict global climate changes responding to variation of cer- 
tain physical parameters |[6T4§l [TTl - n3l[T%ll22l . where the FDT has been used 
largely in its classical formulation 11251 . A vivid demonstration of high pre- 
dictive skill in low-frequency climate response despite structural instability 
of statistical states is given in [12111 . 

Recently, Majda and the author |[3]-[5l developed and tested a novel com- 
putational algorithm for predicting the mean response of nonlinear func- 
tions of states of a chaotic dynamical system to small change in external 
forcing based on the FDT. The major difficulty in this situation is that the 
probability measure in the limit as time approaches infinity in this case is 
typically a Sinai-Ruelle-Bowen probability measure which is supported on 
a large-dimensional (often fractal) set and is usually not absolutely contin- 
uous with respect to the Lebesgue measure [|9j|29l. In the context of Axiom 
A attractors, Ruelle G71l28ll has adapted the classical calculations for FDT 
to this setting. The geometric algorithm (also called the short-time FDT, or 
ST-FDT algorithm in is based on the ideas of OSUSI and takes into 

account the fact that the dynamics of chaotic nonlinear forced-dissipative 
systems often reside on chaotic fractal attractors, where the classical FDT 
formula of the fluctuation-dissipation theorem often fails to produce sat- 
isfactory response prediction, especially in dynamical regimes with weak 
and moderate chaos and slower mixing. It has been discovered in [|3]-[5l 
that the ST-FDT algorithm is an extremely precise response approximation 
for short response times, and can be blended with the classical FDT algo- 
rithm with Gaussian approximation of the state probability density (quasi- 
Gaussian FDT algorithm, or qG-FDT) for longer response times to alleviate 
undesirable effects of expanding Lyapunov directions (which cause numer- 
ical instability in ST-FDT for longer response times). Further developing 
the ST-FDT response algorithm for practical applications, in [|2]| the au- 
thor designed a computationally inexpensive method for ST-FDT using the 
reduced-rank tangent map, and in [1J the ST-FDT algorithm is adapted for 
the response on slow variables of multiscale dynamics, which improves its 
computational stability and simultaneously reduces computational expense. 
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However, dynamical systems describing real-world processes are often 
driven by a stochastic forcing. In this setting, the traditional approach is 
to use the classical FDT algorithm, which computes the linear response to 
small external forcing as a correlation function along a single long-term tra- 
jectory. Typically, it is assumed that the single long-term trajectory samples 
the statistical equilibrium state of the model, however, suitable generaliza- 
tions for dynamics with time-periodic forcing can also be made 11231 12311 . 
A significant drawback of the classical FDT approach is that its compu- 
tational algorithm requires the statistical state probability density together 
with its derivative to be explicitly computed, which is typically not possible 
for complex nonlinear systems. Usually, an approximation is used, such as 
the Gaussian approximation with suitable mean state and covariance ma- 
trix [|3]-[5l. In this case, if the actual statistical state is far from the Gaussian, 
the predicted response is usually considerably different from what is ob- 
served by direct model perturbation (so called ideal response [[3]-[5l). 

On the other hand, the ST-FDT response algorithm is observed to be 
consistently superior to the classical FDT with Gaussian approximation for 
deterministic chaotic dynamical systems with strongly non-Gaussian statis- 
tical states for response times before the numerical instability occurs. In 
this work we adapt the ST-FDT linear response algorithm to be used with 
stochastically forced dynamics (further called stochastic ST-FDT, or SST- 
FDT). Below we observe that the SST-FDT response algorithm, adapted to 
stochastically driven dynamics and blended with the qG-FDT algorithm to 
avoid numerical instability, is also generally superior to the classical FDT 
with Gaussian approximation of the statistical state for both the additive 
and multiplicative noise, just as the ST-FDT algorithm in [3]-[51 for chaotic 
deterministic systems. The manuscript is organized as follows. In Section[2j 
we develop the SST-FDT formula for general time-dependent stochastically 
forced dynamics, and design a practical computational algorithm for au- 
tonomous dynamics with invariant probability measure. In Section[3]we test 
the new algorithm for the stochastically driven Lorenz 96 model [fT9ll20l . 
Section H] summarizes the results of this work. 

2. Fluctuation-dissipation theorem for stochastically driven systems 

Here we consider an Ito stochastic differential equation (SDE) of the form 

(2.1) dx = f a (x, t)dt + o-(jc, t)dW„ 

where x = x(t) e R N , f a : [R N x T] -> R N , a : [R NxK x T] -> R N are 
smooth nonlinear functions, and W t is the ^-dimensional Wiener process. 
Additionally, / depends on a scalar parameter a. We say that the SDE in 
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(12.11) is unperturbed if a = 0, or perturbed otherwise. We also adopt the 
notation / = f Q , with the assumption 

(2.2) ir f a (x,t)\ a= o = B(x)T ] (t), 

oa 

where B(x) is an NxL matrix- valued function, and tj(t) is a L-vector valued 
function. The practical meaning of the above assumption will become clear 
below. 

Let A(jc) be a nonlinear function of x, and let E£„ [A], where t > is the 
elapsed time after to, denote the expectation of A at time to + t over all re- 



alizations of the Wiener process in (12.11) . under the condition that x(to) = x 
(with the short notation Ef [A] = Ef [A]). Let A at the time t be dis- 
tributed according to a probability measure p to , that is, the average value of 
A at time to is 



(2.3) (A)(t ) = Pt0 (A) = f A(x)d Pt0 (x), 

where dp, (x) denotes the measure of the infinitesimal Lebesgue volume dx 
associated with jc. Then, for time to + t, the average of A for the perturbed 
system in (12.11 ) is given by 

(2.4) (A) a (t + t)=p t0 (E ! °i[A])= f Ey a [A]dp t0 (x). 

Jr n ' 

In general, for the same initial distribution p t0 , the average value (A) a (t + 1) 
depends on the value of a. Here, we define the average response 8{A) a (to+t) 
as 

(2.5) 6(A) a (t + t)= f (E^]-E^]Wo(*)- 

The meaning of the average response in (12.51) is the following: for the same 
initial average value of A it provides the difference between the future aver- 
age values of A for the perturbed and unperturbed dynamics in (12.11) . 

If a is small, we can formally linearize (12.51) with respect to a by expand- 
ing in Taylor series around a = and truncating to the first order, obtaining 
the following general linear fluctuation-response formula: 

(2.6) 6(A) a (t + t) = a [ d a Ef [A] dp t0 (x), 

JR N 

where we use the short notation 

(2.7) d a * = ^ 

oa 



a=0 
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2.1. Stochastic short- time linear response. To compute the general lin- 
ear fluctuation-response formula in (12.61 ). we need a suitable algorithm for 
daEY [A]. Let x(t + t) = cpf'x be the trajectory of (ED) starting at x at 
to for a particular realization of the Wiener process W[ t0 ... to +t]- Then, the 
expectation [A] is given by 

(2.8) E*£[A] =E[A(0*''x)], 

where the expectation in the right-hand side is taken with respect to all 
Wiener paths. Therefore, 

(2.9) <9JEf [A] = E [DA (<^ f jt) d a ^x] , 

where DA denotes the derivative of A with respect to its argument. For 
d a (p to,t x, by taking the difference between the perturbed and unperturbed 



versions of (12.11) and linearizing with respect to a at a = 0, we have 

(2 10) dd ^' x = (°f ^°' fjc ' to + t ) dt + D(T (^°' f *> ^ + dW f0+ ,)x 

xd a 4>'°-'x + d a f (<f> toJ x, t Q + t) dt, 

where Df and Der are Jacobians of / and cr, respectively. The above equa- 
tion is a linear stochastic differential equation for d a (f> to,t x with zero initial 
condition (as at to both perturbed and unperturbed solutions start with the 
same jc). It can be solved as follows: let us first introduce the integrating 
factor T'x' (an N x N matrix) given by the solution of the equation 

m = (d/^x, t + t)dt+ 

(2.11) V ' 

+ do- ((p^x, to + 1) dw t0+t )rf, rf = i, 

and represent d a 0' o,t x as a product 

(2.12) d a (f> t 0' t x = T*' t y y, 

where y'°'' is an A/-vector. Then, for the Ito differential of d a (f> toJ x we obtain 
dd a <p toJ x = dTYyY + T'fdyY = (Df((f> to t x, t + t)dt+ 
+Dcr (fx, to + dW t0+t )Ty y y + TfdyY = 
= (Df {4> t0 ' f x, to + t)dt + Da ((f>'°''x, to + t) dW, 0+f )x 

xd a (p ,lh, x + Tf'dyy. 

Comparing the right-hand sides of (12.101) and (12.131) we find that yY satis- 
fies 

(2.14) dyf = {T'fY'dJ^xJo + t)dt, y f f = 0, 
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with the formal solution 

(2.15) yf= f (l*'T 1 d a f(<ffi' T x,t + T)fr. 

Jo 

Therefore, d a (p to,t x is given by 

(2.16) d a ^x= f TfWfT'dJi^x^+r)^. 

Jo 

At this point, observe that the solution T'°'' of (12.111) can be represented as 
a product 

(2.17) T'f = Ty o T ;'- T Tf\ T <t, 

due to the fact that a solution of (|2.1 II) can be multiplied by an arbitrary 
constant matrix on the right and still remains the solution. Then, (12.161) 
becomes 

(2.18) d a <p^x = f T'p'^dJ (cf> t0 ' T x, t + t) dr. 

For smooth f a and cr in (12.11) , (f> to ''x smoothly depends on jc [17], and the 
integrating factor T'° J is in fact the tangent map for the trajectory (f> to ''x: 

q 

(2.19) jy = —^x. 

ox 

With (127181) . (12791) becomes 

(2.20) d a Ef< [A] = JT E[DA (<fP x) T'p; r d a f (^ T x, t + r) ]dr. 
Recalling (12.21) . we write the above formula as 

(2.21) d a Ey [A] = j E[DA (^x) T'lX'B (^ T x) ]ij(t + r)dr. 
Then, the general linear response formula in (12.61) can be written as 

(2.22) 6(A) a (t + 1) = a f R SST (t , t,r)T}(t + T)dr, 

Jo 

where the linear response operator Rssrih, t, r) is given by 

(2.23) R SST (t , t, r) = E f DA (<p^x) T$£*B (^ T x) dp to (x). 

Jr n v 

Further we refer to (12.231) as the stochastic short-time fluctuation-dissipation 
theorem algorithm, or SST-FDT algorithm. The reason is that in practice the 
computation of the tangent map in (12.1 II) for large t becomes numerically 
unstable because of exponential growth due to positive Lyapunov exponents 
(just as observed in 03-[H for deterministic chaotic dynamics). Note that if 
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the stochastic forcing is removed from (12.11) . the SST-FDT response opera- 
tor becomes the usual ST-FDT from [fl3-[5]|. Apparently, (12.231) requires the 
average with respect to p tQ . If p to is not known explicitly, there are some 
opportunities to replace the p-average with time average, particularly for 
the autonomous dynamics with p tQ being the invariant probability measure, 
and also for non- autonomous dynamical systems with explicit time-periodic 
dependence (as done in H23 11241 for classical FDT response). 

2.2. Classical linear response. The standard way to derive the classical 
linear response formula is through the Fokker-Planck equation (or, as it is 
also called, the forward Kolmogorov equation) for the perturbed system in 
(12.11) by neglecting the terms of higher order than the perturbation, as it is 
done in lT3l-[5l l22[|24[|25L However, for the sake of clarity, here we show the 
derivation of the classical FDT directly from (12.61) . Under the assumption 
of continuity of p tQ with respect to the Lebesgue measure, that is, dp, (x) = 
p tQ (x)dx, where p tQ is the probability density, we can also obtain a formal 
general expression for the classical fluctuation-response formula. Using the 
notations 

d i d d \ 

L FP , a (x, t) = -— • (f a (x, f)«) + — 9 — • (o-o- T (x, f)*), 
^2 24, ox \ ax ax / 



£x, a {x, t , t) = Texp l| dr L FP>a (x, t + r) j , 

which are, respectively, the Fokker-Planck and forward Kolmogorov oper- 
ators, we write the expectation Ej^ [A] in the form 



(2.25) 



A(y)£ KtQ! (y, t Q , t)5(x - y)dy = 



I -CkAu, t)A(y)6(x - y)dy = L\ a (x, t , t)A(x), 



where 6(x) is the Dirac delta-function, and the adjoint is taken with respect 
to the standard inner product under the integral. Then, the general response 
formula with dp, (jc) = p, (x)djc becomes 

6(A) a (t + t) = a f d a Ey [A] p tQ (x)dx = 

Jr n 

(2.26) =a I d a £l(x,t ,t)A(x)p t0 (x)dx = 

Jr» 

=a I A(x)d a £ K (x,t ,t)p t0 (x)dx. 

Jr n 
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It is not difficult to show that the parametric derivative of an ordered expo- 
nential of a linear operator L a (x, t) is computed as 



(2.27) 



^TexpQ^ drL a (x,r)j = 

= ^ dr Texp ^ J ds L a (x, s) j ^^'^ Texp |^ ds L a (x, s)j. 

As a result, we obtain 

S(A) a (t + t) = a f dr f £\{x, t + r,t- r)x 
Jo Jr n 

(2.28) xA(x)d a L FP (x, t + t)£ k (x, t , T)p tQ (x)dx = 

= a f dr f E* +r ''- r [A]5 a L FP (x,?o + TK +r (jc)djc, 
Jo Jr w 

where p to+T (x) is given by 

(2.29) p t0+T (x) = £ K (x, t , r)p t0 (x). 

Recalling (12.21) . we recover the classical linear fluctuation-response formula 
in the form 

8(A >„(f + t) = a f dr f ~R'° +TJ ~ T [A] x 

(2.30) Jo JrN 

xd a L FP (x, t Q + T)p t0+T (x)dx = a I R da sAk, t, T)ij(t Q + r)dr, 

Jo 

where the classical linear response operator R c i ass is given by 

(2.31) Rdassito, t, r) = -E f A(^ 0+r - r - r jc)|- • (B(x)p t0+T (x))dx. 

J r n dx 

Observe that, unlike (|2.23l) . in (12.311) one has to know p tQ+T (x) for all re- 
sponse times explicitly to perform differentiation with respect to jc. Usually, 
an approximation is used, such as the Gaussian approximation j3}Q. 

2.3. Special case for autonomous dynamics with ergodic invariant prob- 
ability measure. Here we consider the case where / and er in (12.11) do not 
explicitly depend on t (although f a does with a ^ 0), and we choose p tQ = p 
to be an ergodic invariant probability measure for (12.11) . In this situation, 
one can replace the averaging with respect to the measure p with averaging 
over a single long-term trajectory which starts with an initial condition x in 
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the support of p: 

R SST (t ,t,T) = Elim- f DA(<j) t{ht (p t0 -^x)x 
(2.32) r Jo 

where, without loss of generality, the starting is time t - s, that is, the 
averaging occurs over the endpoints of (f>'° ~ s,s x. Combining the solution 
operators, we obtain 



(2.33) 



R SST (t ,t,r) = Elim- f DA ((f> t0 ~ s,s+t x) x 

r->oc r Jo 

xT^B^-^ds. 



Since the averaging over all independent realizations of the Wiener process 
is needed, we can average over many statistically independent chunks of the 
Wiener path along a single long-time trajectory by setting t = s - r: 

(2.34) R SST (t, r) = lim - f DA (cf>' T ' s+t x) T S 'Z T B (<f>~ T ' s+T x) ds. 

r^oo r J * 

Finally, replacing x with (f>°~ T x (which for finite r is also in the support of 
p), we find that 

(2.35) R SST (t,r) = lim- f DA(/' m - r x)r;^(/'' v x)d^, 
or, denoting jt(s) = o,i x, 

(2.36) R SST (t,T) = Km- f DA(x(s + t - t))T°£7B (x(s))ds. 

r->°° r Jo 

Now, the linear response formula in (12.221) and the response operator in 
(12.231 ) become, respectively, 

6(A) a (t + t) = a f R SST (t - T)T](t + r)dr, 

(2.37) Jo 

tf ssr (0 = lim - DA(x( S + t))r x 'B(x( S ))ds. 

r^oc r Jo xw 

In a similar fashion, for the classical linear response in (12.301) we note that 
the Fokker-Planck operator L FP does not depend on t, and both the forward 
Kolmogorov operator £ K and its adjoint do not depend on to. Taking into 
account that p tQ+T (x) = p(x), where p(x) is the invariant probability density, 
we write 

A^'x)— ■ (B(x)p(x))dx, 

N OX 
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or, after replacing the p-average with the average over the long-term trajec- 
tory, 



Here the expectation can be removed since the averaging over different 
Wiener paths will automatically occur as the long time average is computed. 
As a result, we obtain 



3. Application for the stochastically driven Lorenz 96 model 

The 40-mode deterministic Lorenz 96 model (L96) has been introduced 
by Lorenz and Emanuel HI 9 [1201 as a simple model with large scale features 
of complex nonlinear geophysical systems. The deterministic Lorenz 96 
(L96) model is given by 



with periodic boundary conditions given by X n±N = X n , where N = 40, 
and F being a constant forcing parameter. The model in (13.11 ) is designed 
to mimic midlatitude weather and climate behavior (in particular Rossby 
waves), so periodic boundary conditions are appropriate. It is demonstrated 
in Chapter 2 of [22] that the dynamical regime of the L96 model varies 
with changing the value of constant forcing F: weakly chaotic dynamical 
regimes with F = 5,6, strongly chaotic regime with F = 8, and turbulent 
regimes F = 12, 16, 24 with self-similar time autocorrelation decay. 
Here we apply the stochastic forcing to the L96 model as 



(3.2) dX k = [X k ^(X M - X k _ 2 ) - X k + F] dt + (<r(X)) k (dW,) k , 



where cr : R — > R is a vector-valued function of X, W is a Af-dimensional 
Wiener process, and (dW t ) k is the k-th component of dW (that is, effectively 
cr is a diagonal matrix multiplying the vector dW). As the stochastic Lorenz 
96 (SL96) model above does not depend explicitly on time (except for the 
Wiener noise), we can assume that it has an invariant probability measure 
P- 

In this work, we perturb the SL96 model in (13.21) by a small parameter a 

as 



(3.3) dX k = [X k ^(X k+1 - X k _ 2 ) - X k + F + arj k ] dt + {<r(X)) k {dW t )u, 



(2.39) 




(2.40) 




(3.1) 



X n = Z„_i(Z„ +1 - X„_ 2 ) - X n + F, 1 < k < N, 



IMPROVED LINEAR RESPONSE FOR STOCHASTICALLY DRIVEN SYSTEMS 



11 



where 77 6 R is a constant forcing vector perturbation, which is "turned on" 
at time t = 0. With the invariant probability state p, and the perturbation 
given in (13.31) . the general response formula in (12.61) becomes 



where subscripts for R and 7? are omitted as both the SST-FDT and clas- 
sical response operators apply. We also set the observable A(x) = x, that 
is, the response of the mean state is computed. As an approximation for 
the invariant probability density for the classical response, we choose the 
Gaussian distribution with the same mean and covariance as the actual in- 
variant probability measure, which are determined by averaging along the 
long-term time series of unperturbed (13.21) , and, thus, further call it quasi- 
Gaussian FDT (qG-FDT) as in [3]-[5]|. In this setting, the short-time and 
quasi-Gaussian linear response operators become 



where x and C are the mean state and covariance matrix of the long-time 
series of unperturbed (13.21) . 

3.1. Blended SST/qG-FDT response. Following mm, we also compute 
the blended SST/qG-FDT response as 

(3.6) R SS T/ qG (t) = [1-H(t- ? cutoff )] R S st(!) + H(t- t cut0 $) R qG {t), 

where the blending function H is the Heaviside step-function. The cut-off 
time fcutoff is chosen as 

3 

(3.7) r cu toff = — , 

M 

where A\ is the largest Lyapunov exponent (for details see |[3l|5]|). This 
cut-off time allows to switch to the R qG just before the numerical instability 
occurs in Rsst, and, thus avoid the numerical instability. For constant exter- 
nal forcing and the Heaviside blending step-function the blended response 
operators become 



8{A) a {t) = aK(t)j], 



(3.4) 




(3.5) 




(3.8) 




'cutoff 
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3.2. Computational experiments. Below we perform computational ex- 
periments in the following setting: 

• The number of variables (model size) N = 40 

• Constant forcing F = 6. The L96 model is observed to be weakly 
chaotic in this regime (|3l-l5ll22l. and we would like to compare 
the responses for weakly chaotic deterministic dynamics and the 
stochastically driven dynamics 

• The tangent map T'° J in (12.111 ) is computed in the same fashion as 
in [[Ml 

• Forward Euler numerical scheme with time step At = 0.001 for both 
flXZD) and flO) 

• The linear response is tested for the following settings of the sto- 
chastic term cr: 

- crk = (fully deterministic regime without stochastic forcing) 

- crfc = 1 (additive noise) 

- <Jk = 0.2Xt, cr k = 0.5Xk (multiplicative noise) 

• We compute the linear response operators JZsst, Kqc an d 'Rsst/qg, 
which are given by (13.41) and (13.81) . and compare them with the 
ideal response operator ( Rid ea u which is computed through the direct 
model perturbations Q3-E1 

• The time-averaging is done along a time series of 10000 time units 

• The ideal response operator %deai is computed via direct perturba- 
tions a 10000-member statistical ensemble 

• The comparison of the FDT response operators with the ideal re- 
sponse operator is carried out by evaluating the L 2 relative error 

/0 nN T W-FDT - ft-idealW 

(3.9) L 2 -error = — , 

WKidealW 

and the correlation function 



(3 - 10) Con= \m \wr ir 

where (•, •) denotes the standard Euclidean inner product. Observe 
that the L 2 error shows the general difference between the FDT and 
ideal responses, while the correlation function shows the extent to 
which the responses are collinear (that is, how well the location of 
the response is determined, without considering its magnitude) 

In Figure \T\ we display the L 2 relative errors between the ideal response 
operator and the FDT response operators, together with the intrinsic error 
in the ideal response operator (which is the result of slight nonlinearity in 
the ideal response due to small but finite perturbations). Observe that in 
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Errors for SL96 model, N=40, F=6, a k =0 
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Figure 1. L 2 -errors of the response operators for SL96 
model, N = 40, F = 6. Straight dotted vertical line denotes 
the blending cut-off time for SST/qG-FDT. "Rideai denotes the 
intrinsic error in the ideal response due to slight nonlinearity. 



the fully deterministic regime (F = 6, cr k = 0) the SST-FDT response pro- 
vides a very precise prediction until the time t x 3, and then the errors in 
the SST-FDT grow exponentially rapidly, which is due to the positive Lya- 
punov exponents and numerical instability in the tangent map. On the other 
hand, the qG-FDT response is not precise (reaching about 80% by the time 
t = 1.5), due to the fact that the invariant probability measure associated 
with the deterministic regime is highly non-Gaussian, and most probably 
not continuous with respect to the Lebesgue measure (that is, it does not 
even possess a density). Remarkably, if we look at the stochastically driven 
regimes cr k = 1 (additive noise) and cr k = 0.2X k , cr k = 0.5X k (multiplica- 
tive noise), we see that the behavior of both the SST-FDT and qG-FDT re- 
sponses is qualitatively the same as in the fully deterministic regime, even 
though the dynamics is qualitatively different. Apparently, the level of noise 
in the two stochastically driven regimes cr k = 1 and cr k = 0.2X k is insuf- 
ficient to "smooth out" the invariant probability measure enough for it to 
resemble the Gaussian state and to destabilize the computation of the tan- 
gent map. However, in the cr k = 0.5X k multiplicative noise regime, the 
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Correlations for SL96 model, N=40, F=6, a k =0 Correlations for SL96 model, N=40, F=6, a k =1 



Time 




Time 
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40, F=6, o k =0.5X k 





Figure 2. Correlations of the FDT response operators with 
the ideal response operator for SL96 model, N = 40, F = 6. 
Straight dotted vertical line denotes the blending cut-off time 
for SST/qG-FDT. 



errors in the initial qG-FDT response are reduced to about 40%, which is 
due to the fact that in this regime the invariant probability measure is closer 
to the Gaussian state because of strong noise. The blended SST/qG-FDT 
response yields the lowest errors in all cases, due to its explicit design to 
avoid numerical instability in the SST-FDT algorithm. 

In Figure [2] we show the correlation functions for the same simulations. 
Observe that, although significant L 2 -errors were observed for the qG-FDT 
algorithm for the fully deterministic regime <Tk = 0, its correlations with the 
ideal response are generally on the level of around 0.7, which is remarkable. 
Also, the correlations of the SST-FDT response with the ideal response are 
roughly 1 (nearly perfect correlation) before the numerical instability mani- 
fests itself. As for the blended SST/qG-FDT response, the best correlations 
are achieved in the stochastically forced regimes o> = 1 (additive noise) 
and <Tk = 0.2Xk, cru = 0.5X/, (multiplicative noise), were the correlations do 
not become lower than 0.95 for all response times. For the fully determin- 
istic case 0-^ = the correlations of the blended SST/qG-FDT response are 
about 0.8. 
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Figure 3. Snapshots of the response operators for SL96 
model at T = 1, N = 40, F = 6. 



In addition to displaying the errors and correlations between the FDT re- 
sponse operators and the ideal response operator, in Figures [3]-[5] we show 
the instantaneous snapshots of the linear response operators at times T = 1, 
T = 2 (which are before the SST/qG-FDT cutoff time) and T = 5 (which is 
after the SST/qG-FDT cutoff time). Although the linear response operator 
at a given time is an 40x40 matrix, it has the property of translational invari- 
ance (just like the L96 model itself), and, thus, can be averaged along the 
main diagonal with wrap-around aliasing of rows (or columns) into a single 
vector. These averaged vectors are displayed in Figures [3H21 Observe that 
for the early times of the response 7 = 1,2 the SST/qG-FDT response is 
virtually indistinguishable from the ideal response. As for the qG-FDT re- 
sponse, its best performance is observed in the case of strong multiplicative 
noise cr k = 0.5X/,, where the discrepancies between the qG-FDT and ideal 
response are not much larger than those between the SST/qG-FDT response 
and the ideal response. This is probably the consequence of the fact that the 
strong multiplicative noise changes the invariant probability density of the 
SL96 model to the point where it is relatively close to the Gaussian. For 
other regimes, by the response time T = 2 significant errors develop in the 
qG-FDT response to the right of the main response diagonal. For the longer 
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Figure 4. Snapshots of the response operators for SL96 
model at T = 2, N = 40, F = 6. 



response time T = 5 and all regimes the blended SST/qG-FDT response 
is very similar to the ideal response, while the qG-FDT response again de- 
velops large discrepancies to the right of the main response diagonal for 
<Jk = 0, l,02X k . For the strong multiplicative noise regime, cr fc = 0.5X k , 
and response time T = 5, the qG-FDT yields lower errors than in the other 
regimes, but is still less precise than the SST/qG-FDT response. 



4. Summary 

The classical fluctuation-dissipation theorem, by its design, is suitable 
for computing the linear response for stochastically driven systems, as it as- 
sumes the continuity of the probability measure of the statistical ensemble 
distribution with respect to the Lebesgue measure (which is guaranteed in 
many stochastically driven systems). However, the drawback of the classi- 
cal fluctuation-response formula is that it requires the probability density to- 
gether with its derivative (or their suitable approximations) explicitly in the 
response formula. Unfortunately, for complex systems with many variables 
such an approximation might not be necessarily available with required pre- 
cision. 
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Figure 5. Snapshots of the response operators for SL96 
model at T = 5, N = 40, F = 6. 



In this work, we develop the stochastic short-time fluctuation-dissipation 
formula (SST-FDT) for stochastically driven systems which does not re- 
quire the probability measure of the statistical state of the system to be 
known explicitly. This formula is the analog of the general linear response 
formula Il3l-[5il9ll28ll for chaotic (but not stochastically driven) nonlinear sys- 
tems. We demonstrate that, before the numerical instability due to posi- 
tive Lyapunov exponents occurs, the SST-FDT for the stochastically driven 
Lorenz 96 model is generally superior to the classical FDT formula where 
the probability density of the statistical state is approximated by the Gauss- 
ian density with the same mean and covariance (qG-FDT). We test the new 
SST-FDT formula for the L96 model with stochastic forcing for both the 
additive and multiplicative noise, and observe that the SST-FDT response 
formula is generally better than the qG-FDT in both the error and correla- 
tion comparison, before the numerical instability develops in the SST-FDT 
response. Additionally, the blended SST/qG-FDT response with a simple 
Heaviside blending function clearly performs on top of both the qG-FDT 
and SST-FDT in all studied regimes. The results of this work suggest that 
the SST/qG-FDT algorithm can be used in practical applications with sto- 
chastic parameterization, such as the climate change prediction. 
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